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Abstract 

In  this  paper  we  investigate  the  possibility  of  accelerating  the  transient  simulation  of  MOS 
devices  by  using  waveform  relaxation.  Standard  spatial  discretization  techniques  are 
used  to  generate  a  large,  sparsely-connected  system  of  algebraic  and  ordinary 
differential  equations  in  time.  The  waveform  relaxation  (WR)  algorithm  for  solving  such  a 
system  is  described,  and  several  theoretical  results  that  characterize  the  convergence  of 
WR  for  device  simulation  are  given.  In  addition,  one-dimensional  experimental  results 
are  presented. 
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Abstract 

•' — ^ 

In  this  paper  we  investigate  the  possibility  of  accelerating  the  tran¬ 
sient  simulation  of  MOS  devices  by  using  waveform  relaxation.  Standard 
spatial  discretization  techniques  are  used  to  generate  a  large,  sparsely- 
connected  system  of  algebraic  and  ordinary  differential  equations  in  time. 

The  waveform  relaxation  (WR)  algorithm  for  solving  such  a  system  is  de¬ 
scribed,  and  several  theoretical  results  that  characterize  the  convergence 
of  WR  for  device  simulation  are  given.  In  addition,  one-dimensional  ex¬ 
perimental  results  are  presented. 

r 

1  Introduction 

Both  digital  and  analog  MOS  circuit  designers  rely  heavily  on  circuit  simulation 
programs  like  SPICE  [3]  to  insure  the  correctness  and  to  test  the  performance  of 
their  designs.  For  most  applications,  the  lumped  MOS  models  used  in  these  pro¬ 
grams  [9]  accurately  reflect  the  behavior  of  terminal  currents  and  charges,  but  in 
some  cases,  these  models  are  not  adequate.  In  particular,  charge  redistribution 
between  source  and  drain  during  device  switching  cannot  easily  be  modeled  by 
a  lumped  device,  but  the  details  of  this  charge  redistribution  can  have  an  im¬ 
portant  effect  on  circuit  behavior.  In  circuits  like  dynamic  memory  cells,  sense 
amplifiers,  analog-to-digital  converters,  and  high  frequency  operational  ampli¬ 
fiers.  charge  redistribution  effects  may  not  only  degrade  performance,  but  can 
inhibit  proper  function. 

For  these  critical  applications,  sufficiently  accurate  transient  simulations  can 
be  performed  if,  instead  of  using  a  lumped  model  for  each  transistor,  the  transis- 
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tor  terminal  currents  and  charges  are  computed  by  numerically  solving  the  drift- 
diffusion  based  partial  differential  equation  approximation  for  electron  transport 
in  the  device.  However,  simulating  even  a  few  transistor  circuit  in  this  way  is 
very  computationally  expensive,  because  the  accurate  solution  of  the  transport 
transport  equations  an  MOS  device  requires  a  two  dimensional  mesh  with  more 
than  a  thousand  points. 

In  this  paper  we  investigate  the  possibility  of  accelerating  the  transient  sim¬ 
ulation  of  MOS  devices  by  using  waveform  relaxation.  In  the  next  section  we 
start  by  introducing  the  equations  for  transient  device  simulation.  Then  we  view 
the  result  of  applying  commonly  used  spatial  discretization  techniques  to  these 
equations,  generating  a  large,  sparsely-connected  system  consisting  of  algebraic 
and  ordinary  differential  equations  in  time.  In  Section  3  we  present  the  waveform 
relaxation  algorithm  for  solving  such  a  system,  and  suggest  why  it  may  be  par¬ 
ticularly  efficient.  Several  theoretical  results  that  characterize  the  convergence 
of  the  method  are  presented  in  Section  4,  and  one-dimensional  experimental 
results  are  described  in  section  5.  Finally,  conclusions  and  acknowledgements 
are  given  in  section  6. 

2  Classical  Simulation  Equations 

The  terminal  behavior  of  an  MOS  device  is  well  described  by  the  Poisson  equa¬ 
tion  and  the  electron  current-continuity  equation  [5] 

(VH’  +  q(N  -n)  =  0  (1) 

VJn-q^  =  0  (2) 

In  these  equations  V’  is  the  electrostatic  potential,  q  is  the  magnitude  of  elec¬ 
tronic  charge,  n  is  the  electron  concentration,  and  J„  is  the  electron  current 
density  .V  is  the  net  doping  concentration  given  by  N  —  Np  —  N*  where  Np 
and  Na  are  the  donor  and  acceptor  concentrations. 

The  electron  current  density  is  commonly  approximated  by  the  drift-diffusion 
equation: 

Jn  = -q((in  nVtl' -  DnVn)  (3) 

where  is  the  electron  mobility,  and  Dn  is  the  diffusion  coefficient.  An  equa¬ 
tion  system  with  only  n  and  xi’  as  unknowns  is  derived  by  using  (3)  to  eliminate 
J„  from  (2). 

There  are  a  variety  of  ways  to  spatially  discretize  the  system  of  two  equa¬ 
tions  in  the  two  unknowns  n  and  v.  Given  a  rectangular  two  dimensional  mesh, 
a  common  approach  is  to  use  a  finite-difference  formula  for  the  Poisson  equa¬ 
tion,  and  an  exponentially-fit  finite-difference  formula  for  the  current-continuity 
equation.  For  notational  simplicity,  we  will  assume  that  the  mesh  points  are 
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evenly  spaced  a  distance  I  apart,  so  that  the  discretized  Poisson  equation  at 
each  mesh  point  i  is: 


f5Z^_^)  +  9,2(Ar*"n<)  =  0  (4) 

i 

where  n,, t/j,,  and  Nt  are  the  electron  concentration,  the  potential,  and  the  net 
doping  concentration  at  mesh  point i.  The  summation  is  taken  over  the  nodes 
j  surrounding  t  (four  nodes  for  a  mesh  node  «  not  on  the  boundary,  i.e.  north, 
south,  east,  and  west). 

Under  the  same  assumptions,  and  assuming  constant  mobility,  the  discretized 
current-continuity  equation  with  the  drift-diffusion  approximation  becomes: 

^2  [B(Uj  -  -  B(xn  -  Uj)rii ]  -  ql2  (^nij  =  0  (5) 

where  u,  =  qtli/KT  and  B(x)  =  z/(expx  —  1)  is  the  Bernoulli  function  used 
to  exponentially  fit  the  potential  variation  to  the  electron  concentration  varia¬ 
tion  In  this  equation,  the  Einstein  relation  Dn  =  ( KT/q)pn  has  been  used  to 
eliminate  pn. 

If  there  are  m  mesh  points,  then  the  result  of  applying  the  spatial  discretiza¬ 
tion  to  (1),(2),  and  (3)  is  a  sparse  system  of  m  algebraic  constraints,  represented 
by  (4).  and  a  sparsely  connected  system  of  m  ordinary  differential  equations, 
represented  by  (5). 

3  The  Waveform  Relaxation  Process 

The  standard  approach  used  to  solve  these  two  systems  is  to  discretize  the 
37MO  term  in  (5)  with  a  low  order  integration  method  such  as  backward-Euler 
[ll  The  result  is  a  sequence  of  algebraic  systems  in  2m  unknowns,  each  of  which 
can  be  solved  with  some  variant  of  Newton’s  method  and/or  relaxation.  Another 
approach  is  to  apply  relaxation  directly  to  the  differential  equation  system.  This 
leads  to  a  time  waveform  relaxation  process,  as  given  by  the  following  algorithm 

Although  only  the  Gauss-Jacobi  algorithm  is  presented  for  the  sake  of  no- 
tational  simplicity,  a  Gauss-Seidel  version  could  be  created  by  adjusting  the 
iteration  indexes. 

The  \YR  algorithm  reduces  the  problem  of  simultaneously  solving  m  differ¬ 
ential  equations  and  m  algebraic  equations  to  one  of  iteratively  solving  2m  inde¬ 
pendent  equations.  Each  of  the  m  differential  equations  for  the  n<(t)  waveforms 
can  be  solved  with  a  numerical  integration  method  such  as  backward-Euler. 
Since  they  only  contribute  algebraic  constraints,  the  equations  for  calculating 
the  waveforms  need  to  be  solved  only  at  the  discrete  points  in  time  used 
to  calculate  the  rii(t)  waveforms. 
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Algorithm  1  VVR  Gauss-Jacobi  Algorithm  for  solving  the  system 
produced  by  equations  (4)  and  (5). 

The  superscript  k  denotes  the  iteration  count,  the  subscript 
i  denotes  the  component  index  of  a  vector,  and  e^,  and  e„ 
are  small  positive  numbers. 

k-0 
repeat  { 

k  —k  +  1 

foreach(i  €  { 1 , . . . ,  n } )  { 
solve 

qDn  £  [Biu)-1  -  uf-V}-1  -  -  u-_1) 

-9/2(37"f)  =° 

fot(ik(t).nk(1):  t  £  (0,T].nf(0)  =  ni0) 

} 

}  until(||t’*  -  xl'k  'll  <  e^  and  ||n*  -  ni_1||  <  en) 
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The  inherent  advantage  of  the  WR  approach  is  that  the  differential  equations 
are  solved  in  a  decomposed  fashion,  and  therefore  different  sets  of  timesteps  can 
be  used  at  different  mesh  points  to  calculate  the  time  evolution  of  the  electron 
concentration.  The  method  exploits  multi-rate  behavior.  In  MOS  devices,  the 
rate  at  which  electron  concentrations  evolve  may  be  very  different  in  the  channel 
compared  to  the  source  or  the  drain.  Therefore,  WR  may  prove  to  be  efficient 
for  the  device  simulation  problem,  provided  it  converges,  and  doesn’t  take  too 
many  iterations.  This  is  the  subject  of  the  next  section. 

4  Theoretical  Results 

As  is  usually  the  case  for  waveform  relaxation  algorithms  applied  to  systems  of 
differential  equations.  Algorithm  1  converges  to  the  solution  of  the  differential- 
algebraic  system  for  any  initial  guess  that  matches  the  initial  conditions.  The 
precise  statement  is  given  in  the  following  theorem. 

Theorem  1  Given  a  finite  interval  [0,  T],  and  any  initial  guess  n°(<)  an 
t  €  [0.T],  such  that  n°(0)  =  no,  the  sequence  of  waveforms  produced  by  Alg  1 
converges  to  the  exact  solution  of  the  system  given  by  equations  (4)  and  (5). 

The  proof  of  the  above  theorem  follows  the  same  steps  as  the  Picard-like 
proofs  of  waveform  relaxation  for  ordinary  differential  equations  [10],  First  the 
equations  that  describe  the  difference  between  one  iteration  and  the  next  are 
organized  into  the  form 

6il*i+i  =  Abrfk  +  B6nk(t)  (6) 


*nt+1«)  =  [f(nk+1(t).nk(t)^k(t))  -  /(n*(f),  tf*_1(0)]  (7) 

Jo 

where  Si  k  =  V’f  —  t/’f-1,  Snk  =  nf  —  nf-1.  The  matrices  A,  B  6  3?mxm  and 
the  function  /  :  sftrnxmxm  _  are  constructed  from  the  iteration  equations 
in  Alg  1.  The  next  step  is  to  show  that  (6)  and  (7)  represent  a  contraction. 
To  this  end,  consider  an  interval  of  time  short  enough  to  insure  equation  (7) 
represents  a  contraction  with  respect  to  n  for  a  fixed  t/>.  That  (6)  is  a  contraction 
with  respect  to  1 1-  for  a  fixed  n  is  well-known  [8],  as  (6)  represents  relaxation 
applied  to  the  Poisson  equation.  One  can  fit  the  two  contractions  together  to 
show  that  relaxation  applied  to  the  coupled  system  converges. 

The  above  proof  outline  suggests  that  the  WR  algorithm  converges  in  a 
nonuniform  manner.  That  is.  first  convergence  is  achieved  over  a  short  time 
interval,  set  by  what  is  needed  to  make  (7)  a  contraction,  then  over  the  next 
short  time  interval,  and  then  the  next,  continuing  slowly,  until  the  convergence 
is  achieved  throughout  an  entire  interval  of  interest.  When  applied  to  general 
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differential  equation  systems,  like  circuits,  WR  does  demonstrate  this  nonuni¬ 
formity  in  the  convergence  [7],  but  WR  does  not  usually  show  nonuniformity 
when  applied  to  the  transient  device  simulation  problem. 

In  order  to  analyze  why  this  is  the  case,  we  will  consider  a  model  problem  of 
just  the  differential  equation  associated  with  the  electron  concentration,  n  and 
assume  that  the  potential  ip  is  known.  The  WR  iteration  update  equation  for 
this  case  is  then 

Dn  £  [B(u}  -  u,)n*  -  B(m  -  «>*+1]  -  l2  =  0  (8) 

for  each  i  €  { 1 . .  . .  m} .  Note  that  given  ip,  (8)  is  a  linear  time-varying  differential 
equation  in  n.  For  this  problem  we  have  the  following  theorem: 

Theorem  2  If  at  each  time  t,  ip(t)  is  such  that  the  electric  field  along  any 
vertical  or  horizontal  line  is  either  constant,  or  monotonically  increasing,  then 
(8)  is  a  contraction  in  a  uniform  norm  on  any  finite  interval  [O.T].  That  is, 

™ax[OT]||«n*:+1(0||  <  ■ymaj;[o,T]||6nl(/)||  (9) 

where  ■)  <  1. 

1  he  proof  of  Theorem  2  is  given  in  the  appendix. 

Since  allowing  the  different  differential  equations  to  take  very  different  timesteps 
is  WR's  main  advantage,  if  this  property  were  limited  to  insure  convergence,  the 
WR  algorithm  would  not  be  effective.  Fortunately,  that  the  WR  algorithm  is  a 
contraction  in  a  unform  norm  on  any  interval  implies  that  the  timesteps  used 
to  numerically  integrate  the  differential  equations  are  almost  unconstrained. 
Given  that  the  different  differential  equations  use  different  timesteps,  interpo¬ 
lation  must  be  used  to  communicate  results  between  equations,  and  if  not  done 
carefully  this  can  cause  nonconvergence.  Linear  interpolation  is  certain  not 
cause  problems,  and  therefore  we  have  the  following  theorem  [7]: 

Theorem  3  Let  each  of  the  m  independent  WR  iteration  update  equations 
given  in  (8)  be  solved  numerically  with  backward- Euler,  with  m  different  sets  of 
timesteps.  In  addition,  assume  that  linear  interpolation  is  used  to  derive  values 
for  the  n's  between  ttme  discretization  points.  Then  this  muhiratc  discretized 
WR  algorithm  for  (8)  converges,  regardless  of  the  timestep  selections. 

5  One  Dimensional  Experiments 

Except  for  Theorem  1,  the  above  theoretical  results  only  apply  under  certain 
conditions,  and  are  only  an  indication  that  the  WR  algorithm  may  be  effective. 

In  order  to  verify  that  the  theoretical  results  apply  in  actual  simulation,  a  one¬ 
dimensional  transient  device  simulation  program  was  written  and  applied  to  a 


one-dimensional  approximation  of  an  MOS  device  with  a  conducting  channel. 
The  doping  distribution  for  the  one-dimensional  device  is  given  in  Fig.  1,  where 
the  tick  marks  denote  the  mesh  points.  Potential  and  electron  concentration 
boundary  conditions  were  given  at  x  =  0.0  and  *  =  3.0/i.  The  boundary  values 
for  the  electron  concentration  were  computed  assuming  charge  neutrality  at  the 
“contacts”. 

The  relaxation  process  was  tested  by  first  solving  the  static  problem  with 
zero  volts  across  the  “device”,  and  then  making  a  step  change  of  five  volts. 
Even  with  this  simple  example,  the  variable-by-variable  WR  algorithm  as  given 
in  Alg.  1  was  ineffective.  The  iterates  did  not  converge  in  a  uniform  manner, 
and  they  converged  very  slowly. 

In  order  to  improve  convergence,  rather  than  using  variable-by-variable  de¬ 
composition.  we  partitioned  the  problem  into  blocks  based  on  two  techniques. 
First,  we  associated  the  electron  concentration  at  node  i,  rii(t)  with  the  potential 
Vi(t)  at  that  node.  Then,  in  order  to  try  to  satisfy  the  assumptions  of  Theorem 
2.  we  placed  together  neighboring  nodes  where  we  expected  rapid  changes  in 
the  electric  field.  The  resulting  partitioning  of  the  nodes  are  boxed  in  Fig.  1. 

The  resulting  waveform  iterations  for  the  slowest  converging  variable,  the 
electron  concentration  for  the  mesh  point  where  the  doping  changes  abruptly,  is 
plotted  in  Fig.  2.  As  the  figure  indicates,  with  the  partitioning  just  described, 
the  WR  process  converges  in  just  a  few  iterations  and  the  contraction  is  uniform 
through  time  as  predicted  by  Theorem  2.  The  simulation  was  rerun  with  very 
coarse  timesteps  to  see  the  effects  on  convergence,  and  the  WR  iterations  for  the 
same  node  is  plotted  in  Fig.  3.  As  the  figure  indicates,  using  coarse  timesteps 
does  not  effect  the  overall  convergence,  although  the  convergence  for  small  i  is 
slowed. 

6  Conclusions  and  Acknowledgements 

In  this  paper  we  presented  some  preliminary  results  that  indicate  the  WR  al¬ 
gorithm  may  indeed  be  efficient  for  device  transient  simulation.  In  particular, 
it  was  shown  that  under  conditions  that  can  be  arranged  for  in  practice,  the 
WR  algorithm  is  a  contraction  in  a  uniform  norm  on  any  interval  [0,7].  Also, 
given  these  same  conditions,  the  relaxation  process  will  still  converge  even  if 
very  different  sets  of  timesteps  are  used  for  the  individual  iteration  equations. 
Finally,  we  verified  the  theoretical  results  on  a  one  dimensional  example. 

There  are  several  aspects  of  WR  that  need  to  be  addressed  if  this  method 
it  to  be  efficent  for  two-dimensional  MOS  transient  device  simulation.  Most 
important,  a  general  algorithm  for  blocking  the  device  must  be  developed.  An 
efficent  approach  for  determining  what  discretization  points  to  use  for  the  alge¬ 
braic  constraints  must  be  considered.  In  addition,  the  efficiency  of  WR  methods 
can  also  be  improved  by  refining  the  timesteps  with  iterations,  or  using  a  single 
waveform-Newton  iteration  to  solve  the  nonlinear  WR  iteration  equations. 
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A  Proof  of  Theorem  2 


The  WR  iteration  equations  applied  to  the  model  problem  (8)  can  be  described 
as 

ni+1(<)  =  £>(<)ni+1(<)  +  M(t)nk(t)  (10) 

where  D(t),M(t)  €  &nxn,  and  D(t)  is  negative  diagonal  matrix.  The  assump¬ 
tions  about  the  electric  field  result  in  values  for  the  Bernoulli  functions  such 
that  D(t)  and  M[i)  will  satisfy  the  relation 


!M»(<)||  >  ||m,i(t)||.  (11) 

where  c,  >  0  and  is  strictly  greater  than  zero  for  those  Vs  corresponding  to  the 
mesh  points  next  to  the  boundaries.  Note  that  this  implies 

||D«)-1A/(<)||  <  7  (12) 

for  7  <  1,  for  some  norm  on  2ftnxn  and  for  all  t. 

Given  the  relationship  between  D(t)  and  M(t),  the  WR  algorithm  applied 
to  a  system  of  the  form  of  (13)  will  contract  in  a  uniform  norm.  This  has  been 
shown  for  the  case  when  D(t)  and  M(t)  are  independent  of  t,  using  Laplace 
transforms  [2].  In  the  time  dependent  case,  the  result  can  be  shown  by  examining 
the  difference  between  iteration  k  and  k  ■+■  1  of  (13)  to  get 

6nf+1(t)  =  d,t(<)6nf+1(0  +  ^m,J(<)6n;i(t)  (13) 

i*i 

for  each  mesh  point  i,  where  (t)  =  n*(t)  —  By  assumption.  d„(f)  <  0 

and  ifinffO)  =  0.  Therefore, 

mai[0.T]|<5n{"M(f)|  <  ^  maZ[o,r]l^4y-|ma^[o,T]l^ti*(t)|.  (14) 

Equation  ( 14)  follows  from  the  fact  that  for  all  values  of  6nk+l(t)  on  the  bound¬ 
ary  of  (or  outside)  the  bounded  region  6nf+1(<)  points  back  into  the  bounded 
region  [6], 

Assembling  the  equation  system  from  (14)  results  in 

max[o.T]l<<;nt+1(<)|  <  mar|0,T]|Z?(f)_1  Af(t)|mox[0,T]|6nl(t)|.  (15) 

Then  in  the  norm  for  which  ||D(f)-1  A/(f)||  <  7  <  1.0, 

max[0iT]||«nt+1(<)||  <  7mar[0  r]||(5ni(t)||.  (16) 


which  proves  the  theorem 


